Doubly Truncated Weibull Minimum Distribution (truncweibull_min)

Doubly Truncated Weibull Minimum Distribution (truncweibull_min)#

scipy.stats.truncweibull_min is the Weibull minimum distribution restricted to an interval and re-normalized. Truncation is common when you only observe values within a window (instrument limits, study design, reporting rules) but the underlying phenomenon is still reasonably Weibull.

Learning goals#

  • Write down the PDF/CDF/PPF and understand how truncation re-normalizes a base Weibull.

  • Compute moments using incomplete gamma functions and verify them numerically.

  • Interpret how the parameters (c, a, b, loc, scale) change the shape and the support.

  • Sample from the distribution with a NumPy-only inverse-transform algorithm.

  • Use scipy.stats.truncweibull_min for evaluation, simulation, and fitting.

Notebook roadmap#

  1. Title & classification

  2. Intuition & motivation

  3. Formal definition (PDF/CDF/PPF)

  4. Moments & properties

  5. Parameter interpretation

  6. Derivations (expectation, variance, likelihood)

  7. Sampling & simulation (NumPy-only)

  8. Visualization (PDF, CDF, Monte Carlo)

  9. SciPy integration (scipy.stats.truncweibull_min)

  10. Statistical use cases

  11. Pitfalls

  12. Summary

Prerequisites#

  • comfort with PDF/CDF and expectation

  • basic calculus substitutions

  • familiarity with the Gamma / incomplete Gamma functions (helpful but not required)

import numpy as np

import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import integrate, optimize, special, stats
from scipy.stats import truncweibull_min

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)

# Record versions for reproducibility.
VERSIONS = {"numpy": np.__version__, "scipy": scipy.__version__, "plotly": plotly.__version__}
VERSIONS

{'numpy': '1.26.2', 'scipy': '1.15.0', 'plotly': '6.5.2'}

1) Title & Classification#

Item

Value

Name

truncweibull_min (SciPy: scipy.stats.truncweibull_min)

Type

Continuous

Base distribution

Weibull minimum (weibull_min)

Support (standardized)

\(x \in (a, b]\)

Parameter space

\(c>0\), \(0\le a<b\le\infty\)

Location/scale

loc ∈ ℝ, scale > 0

SciPy defines the truncation bounds in standardized form:

\[ a = \frac{u_\ell - \mathrm{loc}}{\mathrm{scale}}, \qquad b = \frac{u_r - \mathrm{loc}}{\mathrm{scale}}, \]

so the support in the original units is

\[ x \in (u_\ell, u_r] = (\mathrm{loc} + a\,\mathrm{scale},\; \mathrm{loc} + b\,\mathrm{scale}]. \]

Unless stated otherwise, we work in the standardized form (loc=0, scale=1).

2) Intuition & Motivation#

What it models#

The Weibull minimum distribution is a staple model for time-to-failure / lifetime phenomena. A truncated Weibull is appropriate when the physical process is still Weibull-like, but you only observe outcomes inside a window:

  • Left truncation (e.g., detection limits): values below a threshold are unobservable.

  • Right truncation (e.g., instrument saturation, administrative cutoffs): values above a threshold are not recorded.

  • Double truncation: both happen.

Truncation is not the same as censoring: with truncation, observations outside the window never appear in the dataset.

Typical real-world use cases#

  • Reliability testing with a fixed observation period (components only recorded if they fail between times \(u_\ell\) and \(u_r\)).

  • Survival analysis with delayed entry (left truncation) and administrative endpoints (right truncation).

  • Environmental measurements (wind speeds, rainfall intensities) where sensors have detection/saturation limits.

Relations to other distributions#

  • As \(a\to 0\) and \(b\to\infty\), truncweibull_min reduces to the standard Weibull minimum distribution.

  • When \(c=1\), Weibull becomes Exponential, and truncweibull_min becomes a truncated exponential.

  • When \(c=2\) (with a suitable scale), Weibull relates to the Rayleigh distribution.

  • It is a special case of a generic truncated distribution: conditioning a base distribution to lie in \((a,b]\).

3) Formal Definition#

Base Weibull (standardized)#

For shape parameter \(c>0\), the (standardized) Weibull minimum distribution has

\[ f_W(x;c) = c\,x^{c-1} e^{-x^c},\qquad x>0, \]
\[ F_W(x;c) = 1 - e^{-x^c},\qquad x\ge 0. \]

Its survival function is \(S_W(x;c)=e^{-x^c}\).

Truncated Weibull minimum#

Define the normalization constant

\[ Z(c,a,b) = \mathbb{P}(a < W \le b) = F_W(b;c) - F_W(a;c) = e^{-a^c} - e^{-b^c}. \]

Then the truncated density is

\[ f(x; c,a,b) = \frac{f_W(x;c)}{Z(c,a,b)} = \frac{c\,x^{c-1} e^{-x^c}}{e^{-a^c} - e^{-b^c}}, \qquad a < x \le b. \]

The CDF is

\[\begin{split} F(x;c,a,b)= \begin{cases} 0, & x\le a,\\ \dfrac{F_W(x;c)-F_W(a;c)}{F_W(b;c)-F_W(a;c)} =\dfrac{e^{-a^c}-e^{-x^c}}{e^{-a^c}-e^{-b^c}}, & a < x \le b,\\ 1, & x > b. \end{cases} \end{split}\]

A convenient form for the quantile function (PPF) is obtained by solving for \(e^{-x^c}\):

\[ F^{-1}(q;c,a,b) = \Bigl(-\log\bigl((1-q)e^{-a^c} + q e^{-b^c}\bigr)\Bigr)^{1/c},\qquad 0<q<1. \]

Location/scale#

If \(Y\sim\mathrm{truncweibull\_min}(c,a,b)\) in standardized form, then

\[ X = \mathrm{loc} + \mathrm{scale}\,Y \]

has support \((\mathrm{loc}+a\,\mathrm{scale},\; \mathrm{loc}+b\,\mathrm{scale}]\) and

\[ f_X(x)=\frac{1}{\mathrm{scale}} f_Y\!\left(\frac{x-\mathrm{loc}}{\mathrm{scale}}\right). \]
def _logZ_truncweibull_min(c: float, a: float, b: float) -> float:
    """Stable log normalizer log(exp(-a^c) - exp(-b^c)).

    Notes
    -----
    For b = inf, this is simply -a^c.
    """
    if not (c > 0):
        raise ValueError("c must be > 0")
    if not (a >= 0):
        raise ValueError("a must be >= 0")
    if not (b > a):
        raise ValueError("b must be > a")

    A = a**c
    if np.isinf(b):
        return -A

    B = b**c
    # exp(-A) - exp(-B) = exp(-A) * (1 - exp(-(B-A)))
    return -A + np.log1p(-np.exp(-(B - A)))


def truncweibull_min_logpdf(x: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
    """Log-PDF of the standardized truncweibull_min(c, a, b)."""
    x = np.asarray(x, dtype=float)
    logpdf = np.full_like(x, fill_value=-np.inf)
    mask = (x > a) & (x <= b)
    if not np.any(mask):
        return logpdf

    logZ = _logZ_truncweibull_min(c, a, b)
    xm = x[mask]
    logpdf[mask] = np.log(c) + (c - 1.0) * np.log(xm) - xm**c - logZ
    return logpdf


def truncweibull_min_pdf(x: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
    """PDF of the standardized truncweibull_min(c, a, b)."""
    return np.exp(truncweibull_min_logpdf(x, c, a, b))


def truncweibull_min_cdf(x: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
    """CDF of the standardized truncweibull_min(c, a, b)."""
    x = np.asarray(x, dtype=float)
    out = np.zeros_like(x)

    out[x > b] = 1.0
    mask = (x > a) & (x <= b)
    if not np.any(mask):
        return out

    A = a**c
    xm = x[mask]
    delta_x = xm**c - A

    # Numerator: exp(-a^c) - exp(-x^c) = exp(-a^c) * (1 - exp(-(x^c-a^c)))
    num = -np.expm1(-delta_x)

    if np.isinf(b):
        out[mask] = num
        return out

    delta_b = b**c - A
    den = -np.expm1(-delta_b)
    out[mask] = num / den
    return out


def truncweibull_min_ppf(q: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
    """PPF (quantile function) of the standardized truncweibull_min(c, a, b)."""
    q = np.asarray(q, dtype=float)
    if np.any((q < 0) | (q > 1)):
        raise ValueError("q must be in [0, 1]")

    A = a**c
    sa = np.exp(-A)
    sb = 0.0 if np.isinf(b) else np.exp(-(b**c))

    s = (1.0 - q) * sa + q * sb
    return (-np.log(s)) ** (1.0 / c)


# Sanity check against SciPy
c0, a0, b0 = 2.5, 0.3, 2.0
x = np.linspace(a0 + 1e-6, b0, 500)
q = np.linspace(1e-6, 1 - 1e-6, 500)

pdf_err = np.max(np.abs(truncweibull_min_pdf(x, c0, a0, b0) - truncweibull_min.pdf(x, c0, a0, b0)))
cdf_err = np.max(np.abs(truncweibull_min_cdf(x, c0, a0, b0) - truncweibull_min.cdf(x, c0, a0, b0)))
ppf_err = np.max(np.abs(truncweibull_min_ppf(q, c0, a0, b0) - truncweibull_min.ppf(q, c0, a0, b0)))

pdf_err, cdf_err, ppf_err

(3.3306690738754696e-16, 2.220446049250313e-16, 0.0)

4) Moments & Properties#

Raw moments#

For \(k\ge 0\), the \(k\)-th raw moment exists (and is finite) whenever the support avoids \(0\) or the left tail is integrable. For the standardized distribution with truncation \((a,b]\) and \(a\ge 0\):

\[ \mathbb{E}[X^k] = \frac{\int_a^b x^k\,c x^{c-1} e^{-x^c}\,dx}{e^{-a^c} - e^{-b^c}}. \]

With the substitution \(t=x^c\), this becomes an incomplete gamma integral:

\[ \mathbb{E}[X^k] = \frac{\gamma\!\left(1+\frac{k}{c},\,b^c\right) -\gamma\!\left(1+\frac{k}{c},\,a^c\right)}{e^{-a^c} - e^{-b^c}}, \]

where \(\gamma(s,x)\) is the lower incomplete gamma function.

Mean / variance / skewness / kurtosis#

Let \(m_k = \mathbb{E}[X^k]\). Then

  • Mean: \(\mu=m_1\)

  • Variance: \(\sigma^2 = m_2 - m_1^2\)

  • Skewness: \(\gamma_1 = \mathbb{E}[(X-\mu)^3]/\sigma^3\)

  • Excess kurtosis: \(\gamma_2 = \mathbb{E}[(X-\mu)^4]/\sigma^4 - 3\)

MGF / characteristic function#

The moment generating function (MGF) is

\[ M_X(t)=\mathbb{E}[e^{tX}] = \int_a^b e^{tx} f(x)\,dx. \]

  • If \(b<\infty\), then the support is bounded and the MGF exists for all real \(t\).

  • If \(b=\infty\), the MGF existence matches the base Weibull (e.g. for \(c=1\) it only exists for \(t<1\)).

The characteristic function always exists:

\[ \varphi_X(t)=\mathbb{E}[e^{itX}]. \]

Closed forms are generally not simple; numerical quadrature is typical.

Entropy#

The differential entropy is

\[ H(X) = -\int_a^b f(x)\log f(x)\,dx. \]

SciPy provides truncweibull_min.entropy, and we can also compute it numerically.

def truncweibull_min_raw_moment(k: float, c: float, a: float, b: float) -> float:
    """Raw moment E[X^k] in standardized form.

    Uses the formula with lower incomplete gamma:
        E[X^k] = (γ(1+k/c, b^c) - γ(1+k/c, a^c)) / (exp(-a^c) - exp(-b^c)).
    """
    if k < 0:
        raise ValueError("This helper is written for k>=0.")

    s = 1.0 + k / c
    A = a**c
    B = np.inf if np.isinf(b) else b**c

    gamma_s = special.gamma(s)
    num = gamma_s * (special.gammainc(s, B) - special.gammainc(s, A))
    den = np.exp(-A) - (0.0 if np.isinf(B) else np.exp(-B))
    return float(num / den)


def truncweibull_min_moments(c: float, a: float, b: float) -> dict:
    m1 = truncweibull_min_raw_moment(1, c, a, b)
    m2 = truncweibull_min_raw_moment(2, c, a, b)
    m3 = truncweibull_min_raw_moment(3, c, a, b)
    m4 = truncweibull_min_raw_moment(4, c, a, b)

    var = m2 - m1**2
    mu3 = m3 - 3 * m1 * m2 + 2 * m1**3
    mu4 = m4 - 4 * m1 * m3 + 6 * (m1**2) * m2 - 3 * m1**4

    skew = mu3 / (var ** 1.5)
    kurt_excess = mu4 / (var**2) - 3.0

    return {"mean": m1, "var": var, "skew": skew, "kurt_excess": kurt_excess}


c1, a1, b1 = 1.7, 0.2, 2.0
ours = truncweibull_min_moments(c1, a1, b1)
scipy_mean, scipy_var, scipy_skew, scipy_kurt = truncweibull_min.stats(c1, a1, b1, moments="mvsk")

ours, (scipy_mean, scipy_var, scipy_skew, scipy_kurt)

({'mean': 0.8838055783319058,
  'var': 0.19039626228275253,
  'skew': 0.5000215465747762,
  'kurt_excess': -0.585680805798487},
 (0.8838055783319058,
  0.19039626228275253,
  0.5000215465747815,
  -0.585680805798499))

def mgf_numeric(t: float, c: float, a: float, b: float) -> float:
    """Numerical MGF M(t)=E[e^{tX}] via quadrature."""

    def integrand(x: float) -> float:
        return float(np.exp(t * x) * np.exp(truncweibull_min_logpdf(np.array([x]), c, a, b)[0]))

    val, err = integrate.quad(integrand, a, b, limit=400)
    return float(val)


def cf_numeric(t: float, c: float, a: float, b: float) -> complex:
    """Numerical characteristic function φ(t)=E[e^{itX}] via quadrature."""

    def integrand_cos(x: float) -> float:
        lp = truncweibull_min_logpdf(np.array([x]), c, a, b)[0]
        return float(np.cos(t * x) * np.exp(lp))

    def integrand_sin(x: float) -> float:
        lp = truncweibull_min_logpdf(np.array([x]), c, a, b)[0]
        return float(np.sin(t * x) * np.exp(lp))

    re, _ = integrate.quad(integrand_cos, a, b, limit=400)
    im, _ = integrate.quad(integrand_sin, a, b, limit=400)
    return complex(re, im)


def entropy_numeric(c: float, a: float, b: float) -> float:
    """Numerical differential entropy -E[log f(X)]."""

    def integrand(x: float) -> float:
        lp = truncweibull_min_logpdf(np.array([x]), c, a, b)[0]
        if not np.isfinite(lp):
            return 0.0
        return float(-np.exp(lp) * lp)

    val, _ = integrate.quad(integrand, a, b, limit=400)
    return float(val)


t_values = [-2.0, -1.0, 0.5]
mgf_vals = [mgf_numeric(t, c1, a1, b1) for t in t_values]
cf_vals = [cf_numeric(t, c1, a1, b1) for t in t_values]

H_scipy = float(truncweibull_min.entropy(c1, a1, b1))
H_num = entropy_numeric(c1, a1, b1)

mgf_vals, cf_vals, H_scipy, H_num

([0.23492536535757963, 0.45106303854527163, 1.5944101662659167],
 [(-0.08811862717309919-0.6682569212463516j),
  (0.5811634487536085-0.698157117953174j),
  (0.8829961202882072+0.4168074733874823j)],
 0.4714474041561004,
 0.4714474041561004)

5) Parameter Interpretation#

Shape parameter c (Weibull shape)#

In the (untruncated) Weibull minimum distribution:

  • \(c<1\) produces a decreasing density on \((0,\infty)\) with a heavy right tail (sub-exponential).

  • \(c=1\) reduces to the exponential distribution.

  • \(c>1\) produces a density that rises from 0 and then decays; the hazard rate is increasing.

After truncation, the same tendencies remain, but the distribution is confined to \((a,b]\).

Truncation parameters a and b#

  • Increasing a removes the left part of the support and shifts mass to the right.

  • Decreasing b removes the upper tail and forces the distribution to concentrate below b.

  • Limits:

    • \(b\to\infty\) gives a left-truncated Weibull.

    • \(a\to 0\) and \(b\to\infty\) recovers the base Weibull.

loc and scale#

  • loc shifts the distribution.

  • scale stretches/compresses it.

  • The truncation bounds in the original units are \(u_\ell = \mathrm{loc}+a\,\mathrm{scale}\) and \(u_r = \mathrm{loc}+b\,\mathrm{scale}\).

def plot_pdf_family_by_c(a: float, b: float, c_values: list[float]) -> go.Figure:
    x = np.linspace(a + 1e-6, b, 800)
    fig = go.Figure()
    for c in c_values:
        fig.add_trace(go.Scatter(x=x, y=truncweibull_min_pdf(x, c, a, b), mode="lines", name=f"c={c}"))
    fig.update_layout(
        title="PDF changes with shape c (fixed a,b)",
        xaxis_title="x",
        yaxis_title="density",
        legend_title="shape",
    )
    return fig


def plot_pdf_family_by_trunc(c: float, trunc_pairs: list[tuple[float, float]]) -> go.Figure:
    fig = go.Figure()
    for a, b in trunc_pairs:
        x = np.linspace(a + 1e-6, b, 800)
        fig.add_trace(go.Scatter(x=x, y=truncweibull_min_pdf(x, c, a, b), mode="lines", name=f"a={a}, b={b}"))
    fig.update_layout(
        title="PDF changes with truncation (fixed c)",
        xaxis_title="x",
        yaxis_title="density",
        legend_title="(a,b)",
    )
    return fig


fig1 = plot_pdf_family_by_c(a=0.2, b=2.0, c_values=[0.7, 1.0, 1.7, 3.0])
fig2 = plot_pdf_family_by_trunc(c=1.7, trunc_pairs=[(0.0, 2.0), (0.2, 2.0), (0.2, 1.2), (0.8, 2.0)])
fig1.show()
fig2.show()

6) Derivations#

(a) Raw moment and expectation#

Let \(X\sim\mathrm{truncweibull\_min}(c,a,b)\) in standardized form. For \(k\ge 0\):

\[ \mathbb{E}[X^k] = \frac{1}{Z(c,a,b)}\int_a^b x^k\,c x^{c-1} e^{-x^c}\,dx = \frac{1}{Z(c,a,b)}\int_a^b c x^{k+c-1} e^{-x^c}\,dx. \]

Substitute \(t=x^c\) so that \(dt=c x^{c-1}dx\) and \(x^k=t^{k/c}\):

\[ \int_a^b c x^{k+c-1} e^{-x^c}\,dx = \int_{a^c}^{b^c} t^{k/c} e^{-t}\,dt = \gamma\!\left(1+\frac{k}{c}, b^c\right) - \gamma\!\left(1+\frac{k}{c}, a^c\right). \]

Dividing by \(Z(c,a,b)=e^{-a^c}-e^{-b^c}\) yields the raw moment formula used above.

(b) Variance#

Once \(m_1=\mathbb{E}[X]\) and \(m_2=\mathbb{E}[X^2]\) are available,

\[ \mathrm{Var}(X)=m_2 - m_1^2. \]

(c) Likelihood#

Given i.i.d. observations \(x_1,\dots,x_n \in (a,b]\), the log-likelihood (standardized form) is

\[ \ell(c,a,b\mid x) = \sum_{i=1}^n \log f(x_i;c,a,b) = n\log c + (c-1)\sum_{i=1}^n\log x_i - \sum_{i=1}^n x_i^c - n\log\bigl(e^{-a^c}-e^{-b^c}\bigr). \]

With loc/scale, use \(y_i=(x_i-\mathrm{loc})/\mathrm{scale}\) and add the Jacobian term \(-n\log(\mathrm{scale})\).

def loglike_c_given_ab(c: float, x: np.ndarray, a: float, b: float) -> float:
    """Log-likelihood as a function of c, with (a,b) fixed in standardized form."""
    if c <= 0:
        return -np.inf
    x = np.asarray(x, dtype=float)
    if np.any((x <= a) | (x > b)):
        return -np.inf

    logZ = _logZ_truncweibull_min(c, a, b)
    n = x.size
    return float(n * np.log(c) + (c - 1.0) * np.sum(np.log(x)) - np.sum(x**c) - n * logZ)


# Example: MLE for c when (a,b) are known
true_c, true_a, true_b = 1.8, 0.2, 2.0
data = truncweibull_min.rvs(true_c, true_a, true_b, size=2000, random_state=rng)

def nll(c: float) -> float:
    return -loglike_c_given_ab(c, data, true_a, true_b)

res = optimize.minimize_scalar(nll, bounds=(0.1, 10.0), method="bounded")
c_hat = float(res.x)

c_hat, true_c

(1.7778192822675947, 1.8)

7) Sampling & Simulation#

Because we have a closed-form PPF, sampling is straightforward via inverse transform.

From the truncated CDF, one convenient identity is

\[ e^{-X^c} = (1-U)e^{-a^c} + U e^{-b^c},\qquad U\sim\mathrm{Unif}(0,1). \]

So the algorithm is:

  1. Sample \(U\sim\mathrm{Unif}(0,1)\).

  2. Compute \(S = (1-U)e^{-a^c} + U e^{-b^c}\).

  3. Return \(X = (-\log S)^{1/c}\).

This uses only log, exp, and power, so it is easy to implement with NumPy.

def truncweibull_min_rvs_numpy(
    c: float,
    a: float,
    b: float,
    size: int,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """NumPy-only sampling for standardized truncweibull_min(c, a, b)."""
    if rng is None:
        rng = np.random.default_rng()
    u = rng.random(size)
    return truncweibull_min_ppf(u, c, a, b)


c2, a2, b2 = 1.7, 0.2, 2.0
s_numpy = truncweibull_min_rvs_numpy(c2, a2, b2, size=20000, rng=rng)
s_scipy = truncweibull_min.rvs(c2, a2, b2, size=20000, random_state=rng)

np.mean(s_numpy), np.mean(s_scipy), truncweibull_min.mean(c2, a2, b2)

(0.8850091193506309, 0.8853549491974178, 0.8838055783319058)

8) Visualization#

We’ll visualize:

  • PDF (analytic) and histogram of Monte Carlo samples

  • CDF (analytic) and empirical CDF

c3, a3, b3 = 1.7, 0.2, 2.0
samples = truncweibull_min_rvs_numpy(c3, a3, b3, size=5000, rng=rng)

x = np.linspace(a3 + 1e-6, b3, 600)
pdf = truncweibull_min_pdf(x, c3, a3, b3)
cdf = truncweibull_min_cdf(x, c3, a3, b3)

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF + histogram", "CDF + ECDF"))

# PDF + histogram
fig.add_trace(go.Histogram(x=samples, histnorm="probability density", nbinsx=40, name="samples"), row=1, col=1)
fig.add_trace(go.Scatter(x=x, y=pdf, mode="lines", name="pdf"), row=1, col=1)

# CDF + ECDF
ecdf_x = np.sort(samples)
ecdf_y = np.arange(1, ecdf_x.size + 1) / ecdf_x.size
fig.add_trace(go.Scatter(x=x, y=cdf, mode="lines", name="cdf"), row=1, col=2)
fig.add_trace(go.Scatter(x=ecdf_x, y=ecdf_y, mode="lines", name="ecdf"), row=1, col=2)

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(title=f"truncweibull_min(c={c3}, a={a3}, b={b3})", bargap=0.05)
fig.show()

9) SciPy Integration#

scipy.stats.truncweibull_min provides:

  • pdf, logpdf, cdf, ppf, rvs

  • stats and moment

  • entropy

  • fit (MLE-based parameter estimation)

A common workflow is to freeze the distribution with parameters and then call methods on the frozen object.

c4, a4, b4 = 1.8, 0.2, 2.0
rv = truncweibull_min(c4, a4, b4)  # frozen standardized distribution

x0 = np.array([0.25, 0.5, 1.0, 1.8])
rv.pdf(x0), rv.cdf(x0), rv.ppf([0.1, 0.5, 0.9])

(array([0.597198, 0.847306, 0.72325 , 0.176502]),
 array([0.027815, 0.213994, 0.631769, 0.972303]),
 array([0.357434, 0.830846, 1.51064 ]))

# Fitting: often (a,b) are known measurement bounds, so we fit c with a,b fixed.
sim = rv.rvs(size=3000, random_state=rng)
c_hat_fit, a_hat_fit, b_hat_fit, loc_hat_fit, scale_hat_fit = truncweibull_min.fit(sim, fa=a4, fb=b4, floc=0, fscale=1)
c_hat_fit, (a_hat_fit, b_hat_fit, loc_hat_fit, scale_hat_fit)

(1.787695312500002, (0.2, 2.0, 0, 1))

10) Statistical Use Cases#

(a) Hypothesis testing#

A simple and interpretable test is whether the data are compatible with an exponential tail on \((a,b]\). Because Weibull with \(c=1\) is exponential, we can test

\[ H_0: c = 1 \quad\text{(truncated exponential)} \qquad\text{vs}\qquad H_1: c \ne 1. \]

We can use a likelihood ratio test (LRT) when \((a,b)\) are fixed/known.

(b) Bayesian modeling#

When you want full uncertainty quantification, treat \(c\) as a parameter with a prior (e.g. a Gamma prior) and compute the posterior.

(c) Generative modeling#

Once parameters are estimated, the distribution is a handy bounded generator for synthetic lifetimes or severities (e.g. simulation studies, stress testing, Monte Carlo pipelines).

# (a) Likelihood ratio test for H0: c=1 vs H1: c free (a,b known)
a5, b5 = 0.2, 2.0
x_obs = truncweibull_min.rvs(1.8, a5, b5, size=800, random_state=rng)

c_mle = float(optimize.minimize_scalar(lambda c: -loglike_c_given_ab(c, x_obs, a5, b5), bounds=(0.1, 10.0), method="bounded").x)
ll_alt = loglike_c_given_ab(c_mle, x_obs, a5, b5)
ll_null = loglike_c_given_ab(1.0, x_obs, a5, b5)

lrt = 2.0 * (ll_alt - ll_null)
p_value = float(stats.chi2.sf(lrt, df=1))

dict(c_mle=c_mle, lrt=lrt, p_value=p_value)

{'c_mle': 1.904282724745143,
 'lrt': 83.81779762028248,
 'p_value': 5.425417704887916e-20}

# (b) Simple Bayesian inference for c with a,b fixed: grid posterior
c_grid = np.linspace(0.2, 5.0, 600)

# Prior: Gamma(shape=2, scale=1) on c (mean=2). Feel free to adjust.
log_prior = stats.gamma(a=2.0, scale=1.0).logpdf(c_grid)
log_like = np.array([loglike_c_given_ab(c, x_obs, a5, b5) for c in c_grid])

log_post_unnorm = log_prior + log_like
log_post = log_post_unnorm - special.logsumexp(log_post_unnorm) - np.log(c_grid[1] - c_grid[0])
post = np.exp(log_post)

c_mean = float(np.trapz(c_grid * post, c_grid))

# 95% credible interval by CDF on the grid
cdf_post = np.cumsum(post) * (c_grid[1] - c_grid[0])
c_lo = float(np.interp(0.025, cdf_post, c_grid))
c_hi = float(np.interp(0.975, cdf_post, c_grid))

fig = go.Figure()
fig.add_trace(go.Scatter(x=c_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=c_mean, line_dash="dash", annotation_text=f"mean={c_mean:.2f}")
fig.add_vrect(x0=c_lo, x1=c_hi, opacity=0.15, annotation_text="95% CI")
fig.update_layout(title="Posterior over c (a,b fixed)", xaxis_title="c", yaxis_title="density")
fig.show()

dict(c_mean=c_mean, c_95=(c_lo, c_hi))

{'c_mean': 1.897811196254147, 'c_95': (1.7183616784367284, 2.0639336197638096)}

# (c) Generative modeling: fit c, then generate synthetic samples
c_fit = float(truncweibull_min.fit(x_obs, fa=a5, fb=b5, floc=0, fscale=1)[0])
rv_fit = truncweibull_min(c_fit, a5, b5)

synthetic = rv_fit.rvs(size=5000, random_state=rng)

fig = px.histogram(
    x_obs,
    nbins=40,
    histnorm="probability density",
    opacity=0.55,
    labels={"value": "x"},
    title=f"Observed vs synthetic (fit c={c_fit:.2f}, a={a5}, b={b5})",
)
fig.add_trace(go.Histogram(x=synthetic, nbinsx=40, histnorm="probability density", opacity=0.55, name="synthetic"))

x_line = np.linspace(a5 + 1e-6, b5, 600)
fig.add_trace(go.Scatter(x=x_line, y=rv_fit.pdf(x_line), mode="lines", name="fitted pdf"))
fig.update_layout(barmode="overlay", xaxis_title="x", yaxis_title="density")
fig.show()

11) Pitfalls#

  • Parameter validity: require c>0, a>=0, b>a, and scale>0.

  • Standardized bounds: in SciPy, a and b are standardized; the actual cutoffs are loc + a*scale and loc + b*scale.

  • Data outside support: any observation \(x\le a\) or \(x>b\) has likelihood 0 under the model.

  • Numerical stability: the normalizer \(e^{-a^c}-e^{-b^c}\) can suffer cancellation when \(b^c-a^c\) is tiny; prefer log-space or expm1-based forms.

  • Fitting: if you let fit estimate a and b, it may push them toward sample min/max; in many applications truncation bounds are known and should be fixed.

  • MGF existence when \(b=\infty\): for infinite upper bound, the MGF can diverge for \(t>0\) depending on c (like the base Weibull).

12) Summary#

  • truncweibull_min is a Weibull minimum conditioned to lie in \((a,b]\) (with a,b standardized in SciPy).

  • The PDF is the Weibull PDF divided by \(Z=e^{-a^c}-e^{-b^c}\); the CDF and PPF follow in closed form.

  • Raw moments reduce to incomplete gamma differences, enabling mean/variance/skewness/kurtosis calculations.

  • Sampling is easy via the PPF and can be implemented with NumPy only.

  • For inference, scipy.stats.truncweibull_min supports evaluation, simulation, and MLE-based fitting; consider fixing truncation bounds when they are design constraints.

Contents